#delimit;
clear;
set more off;
version 11;
  quietly log;
  local logon = r(status);
 if "`logon'" == "on" {; log close; };
log using NiemanGibler_JOP_demdiff.log, text replace;

/* Doug Gibler and Mark Nieman */
/* Data Analysis for Taking Democratic Differences Seriously */
/* Input: NiemanGibler_demdiff_ally.dta, NiemanGibler_demdiff_trade */
/* Output:  Table 1, Figure 1, Neglible effect test, Figure 2, Figure 3 */

/* Preamble */
clear;
estimates clear;
set memory 4g;
program drop _all;
clear matrix;
tempfile full;
set seed 17;

	/* Create the split-sample logit */
		#delimit;
		program define ssp_lf, rclass;		
			args lnf beta gamma ;
			tempvar rel violate;
				quietly gen double `rel' = 1/(1+exp(-`gamma'));
				quietly gen double `violate' = 1/(1+exp(-`beta'));
			quietly replace `lnf' = ln((1-`rel')+(`rel'*(1-`violate'))) if $ML_y1==0;	
			quietly replace `lnf' = ln((`rel')*(`violate')) if $ML_y1==1;
		end;
		
	/* Calculate Rubin formula for uncertainty */
		capture program drop rubin;
		program define rubin;
			local names: colnames(Q1);
			/* MI param estimates */
			matrix Q = (Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q9 + Q10)/10;
			/* MI covariances     */
			matrix W = (W1 + W2 + W3 + W4 + W5 + W6 + W7 + W8 + W9 + W10)/10;		
			matrix ll = (ll1 + ll2 + ll3 + ll4 + ll5 + ll6 + ll7 + ll8 + ll9 + ll10)/10;
			local  k = colsof(Q);
			forvalues i=1/10 {;
			  matrix QQ=Q`i'-Q;
			  if `i'==1 {;
				matrix B=(QQ)'*QQ;
				};
			  else matrix B = B + (QQ)'*QQ;
			};
			/* Covariance adjustment */
			matrix B=B/(10-1);		
			/* Final covariance   */
			matrix V=W+(1+1/10)*B;				
			ereturn post Q V ;
		end;

/* Load alliance data */
	use NiemanGibler_demdiff_ally, clear;

/* Table 1 */
	/* Replicate Leeds et al (Model 1) and identify shared observations*/
	#delimit;
	logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
		asymmetry  nomicoop treaty milinst time timesquare timecube
		, robust cluster(atopidphase);
	estimates store LMV;
	gen LMV=1 if e(sample)==1; 
	matrix lmv_b = e(b);
	matrix lmv_se = e(V);

	/* Run Relevancy equation to identify shared observations*/
	logit violate TerrThreat2 num_dems_on_border  num_borders  cold 
		Lenergy_pop Oil_exp  maj_pow kgd_riv , cluster(atopidphase) robust;
	gen GT=1 if e(sample)==1;
	matrix B1 = e(b);

	/* Idenity cases in both samples*/
	gen GL=0;
	replace GL=1 if GT==1 & LMV==1;
	/* Save temp file to use later */
	save `full', replace;

	/* Rerun on Leeds et al on same sample (GL) (Model 2)*/
	logit violate wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
		asymmetry  nomicoop treaty milinst time timesquare timecube if GL==1
		, robust cluster (atopidphase);
	estimates store LMV2;
	matrix B2 = e(b);
	matrix ll_lmv = e(ll);

	/* Model 3*/
	/* Split-sample model applying Rubin for uncertainty in estimates */
	forvalues i=1/10 {;
	ml model lf ssp_lf  
		(B: violate = wcchangedA capcheither20 regchangeeither chTEAltorBmt alformA demdA 
		asymmetry  nomicoop treaty milinst time timesquare timecube)
		(G: violate = TerrThreat`i' num_dems_on_border   num_borders  
		 cold Lenergy_pop Oil_exp maj_pow kgd_riv  ) if GL==1 
		, robust cluster(atopidphase) tech(nr);
	ml init B2 B1 , copy ;
	ml search;
	ml maximize, diff;
	matrix Q`i' = e(b);
		matrix W`i' = e(V);
		matrix  ll`i' = e(ll);
		matrix chi`i' = e(chi2);
		nois _dots `i' 0;
	};

	rubin;
	/* -estout- is necessary to use -estimates store- */	
		estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
		matrix rename ll ll_full;
		matrix fitB = e(b);
		matrix fitSE = e(V);
	estimates store full;	
	
	/* Table 1 results */
	#delimit;
	estout LMV LMV2 full, stats(ll N )
	cells(b(fmt(3) star) se(par))
	modelwidth(8)
	starlevels(* 0.05) legend
	varwidth(15)
	style(tex) replace;	
	matrix list ll_full;
	
/* Figure 1: Density of Territorial Threat and Frequency of Violations.*/		
	keep if GL==1;	
	egen terrthreat_bin = cut(Threat_combine), at(0(0.005)1);
	egen tot_vio = total(violate), by(terrthreat_bin);
	#delimit;
	twoway kdensity Threat_combine || bar tot_vio terrthreat_bin, barw(.007) scheme(s1mono) 
		legend(lab(1 "Density of Observations") lab(2 "Violations") symx(*.5))
		xlabel(0(.1).7)
		ytitle("Number of Violations") xtitle("Territorial Threat") 
		note("Note: The frequency of violations is overlayed with the kernel density of Territorial Threat within the"
			"sample from Table 1 Model 3.");
	graph save fig1, replace asis;

/*Rainey test for neglible effect */
	/* Collapse data to mean/median and expand for simulations */
	collapse  	(mean)  Lenergy_pop time num_dems_on_border num_borders 
				(median)  
					demdA Oil_exp cold  
							nomicoop treaty asymmetry milinst regchangeeither chTEAltorBmt capcheither20
							wcchangedA alformA maj_pow kgd_rival
					Threat_combine;
	expand 10000;
	gen timesquare = time^2; gen timecube = time^3;
	drawnorm beta1 beta2 beta3 beta4 beta5 beta6 beta7 beta8 beta9 beta10 beta11 beta12 beta13 beta14
	gamma1  gamma3 gamma4 gamma5 gamma6 gamma7 gamma8 gamma9 gamma10, 
		means(fitB) cov(fitSE);
	drawnorm beta1_lmv beta2_lmv beta3_lmv beta4_lmv beta5_lmv beta6_lmv beta7_lmv 
		beta8_lmv beta9_lmv beta10_lmv beta11_lmv beta12_lmv beta13_lmv beta14_lmv, 
		means(lmv_b) cov(lmv_se);
	/* Generate latent value for relevancy equation */
	gen double gbhat = gamma1*Threat_combine  + gamma3*num_dems_on_border + gamma4*num_borders + gamma5*cold +
						gamma6*Lenergy_pop + gamma7*Oil + gamma8*maj_pow + gamma9*kgd_rival + gamma10;
	/* Generate latent value for outcome equation, varying change in leader's societal coalition */
	gen double xbhat_0 = beta1*0 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*demdA + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	gen double xbhat_1 = beta1*1 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*demdA + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	/* Generate latent value for replication of Leeds et al, varying change in leader's societal coalition */
	gen double xbhat_lmv_0 = beta1_lmv*0 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*demdA + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;
	gen double xbhat_lmv_1 = beta1_lmv*1 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*demdA + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;
	/* Generate latent value for outcome equation, varying democracy */			
	gen double xbhat_d0 = beta1*0 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*0 + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	gen double xbhat_d1 = beta1*0 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*1 + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	/* Generate latent value for replication of Leeds et al, varying democracy */	
	gen double xbhat_lmv_d0 = beta1_lmv*0 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*0 + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;
	gen double xbhat_lmv_d1 = beta1_lmv*0 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*1 + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;			
	/* Generate predicted probabilities for change in leader's societal coalition 
	for each condition in both the split-sample and replication equations */				
	gen double G = 1/(1+exp(-gbhat));
	gen double F_0 = 1/(1+exp(-xbhat_0));
	gen double F_1 = 1/(1+exp(-xbhat_1));
	gen double vio_0 = G*F_0;
	gen double vio_1 = G*F_1;
	gen double lmv_0 = 1/(1+exp(-xbhat_lmv_0));
	gen double lmv_1 = 1/(1+exp(-xbhat_lmv_1));
	/* Generate predicted probabilities for democracy for each condition in both
	the split-sample and replication equations */	
	gen double vio_d1 = G*(1/(1+exp(-xbhat_d1)));
	gen double vio_d0 = G*(1/(1+exp(-xbhat_d0)));
	gen double lmv_d0 = 1/(1+exp(-xbhat_lmv_d0));
	gen double lmv_d1 = 1/(1+exp(-xbhat_lmv_d1));				
	gen double lmv_fd_d = lmv_d1 - lmv_d0;
	gen double vio_fd_d = vio_d1 - vio_d0;
	gen double lmv_fd_sc = lmv_1 - lmv_0;
	gen double vio_fd_sc = vio_1 - vio_0;
	/* Calculate mean and standard devation of predicted probabilities, removing extreme values */				
	collapse (mean) lmv_fd_d vio_fd_d lmv_fd_sc vio_fd_sc 
			 (p2) lmv_lo_d=lmv_fd_d vio_lo_d=vio_fd_d lmv_lo_sc=lmv_fd_sc vio_lo_sc=vio_fd_sc 
			 (p98) lmv_hi_d=lmv_fd_d vio_hi_d=vio_fd_d lmv_hi_sc=lmv_fd_sc vio_hi_sc=vio_fd_sc;

	sum lmv_lo_d lmv_fd_d lmv_hi_d vio_lo_d vio_fd_d vio_hi_d
		lmv_lo_sc lmv_fd_sc lmv_hi_sc vio_lo_sc vio_fd_sc vio_hi_sc;	

/* Figure 2: Predicted Probabilities of an Alliance Violation, Democracy, and Territorial Threat */
	/* Load data tempfile */
	use `full', clear;	
	/* Collapse data to mean/median and expand for simulations */
	collapse  	(mean)  Lenergy_pop time num_dems_on_border num_borders 
				(median)  
					demdA Oil_exp cold  
							nomicoop treaty asymmetry milinst regchangeeither chTEAltorBmt capcheither20
							wcchangedA alformA maj_pow kgd_rival
				(p2) terrthreat_min=Threat_combine (p98) terrthreat_max=Threat_combine;
		expand 100;
		gen terrthreat = terrthreat_min + (terrthreat_max - terrthreat_min)*(_n-1)/_N;
		expand 100;
		gen timesquare = time^2; gen timecube = time^3;
		drawnorm beta1 beta2 beta3 beta4 beta5 beta6 beta7 beta8 beta9 beta10 beta11 beta12 beta13 beta14
		gamma1  gamma3 gamma4 gamma5 gamma6 gamma7 gamma8 gamma9 gamma10, 
				means(fitB) cov(fitSE);
		drawnorm beta1_lmv beta2_lmv beta3_lmv beta4_lmv beta5_lmv beta6_lmv beta7_lmv beta8_lmv beta9_lmv beta10_lmv beta11_lmv beta12_lmv beta13_lmv beta14_lmv, means(lmv_b) cov(lmv_se);
	/* Generate latent value for relevancy equation */
	gen double gbhat = gamma1*terrthreat  + gamma3*num_dems_on_border + gamma4*num_borders + gamma5*cold +
						gamma6*Lenergy_pop + gamma7*Oil + gamma8*maj_pow + gamma9*kgd_rival + gamma10;
	/* Generate latent value for outcome equation, varying change in leader's societal coalition */
	gen double xbhat_0 = beta1*0 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*demdA + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	gen double xbhat_1 = beta1*1 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*demdA + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	/* Generate latent value for replication of Leeds et al, varying change in leader's societal coalition */
	gen double xbhat_lmv_0 = beta1_lmv*0 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*demdA + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;
	gen double xbhat_lmv_1 = beta1_lmv*1 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*demdA + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;
	/* Generate latent value for outcome equation, varying democracy */		
	gen double xbhat_d0 = beta1*0 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*0 + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	gen double xbhat_d1 = beta1*0 + beta2*capcheither20 + beta3*regchangeeither + beta4*chTEAltorBmt + beta5*alformA + beta6*1 + 
						beta7*asymmetry + beta8*nomicoop + beta9*treaty + beta10*milinst + beta11*time + beta12*timesquare + beta13*timecube + beta14;
	/* Generate latent value for replication of Leeds et al, varying democracy */	
	gen double xbhat_lmv_d0 = beta1_lmv*0 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*0 + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;
	gen double xbhat_lmv_d1 = beta1_lmv*0 + beta2_lmv*capcheither20 + beta3_lmv*regchangeeither + beta4_lmv*chTEAltorBmt + beta5_lmv*alformA + 
					beta6_lmv*1 + beta7_lmv*asymmetry + beta8_lmv*nomicoop + beta9_lmv*treaty + beta10_lmv*milinst + beta11_lmv*time + 
					beta12_lmv*timesquare + beta13_lmv*timecube + beta14_lmv;			
					
	/* Generate predicted probabilities for change in leader's societal coalition 
	for each condition in both the split-sample and replication equations */
	gen double G = 1/(1+exp(-gbhat));
	gen double vio_0 = G*(1/(1+exp(-xbhat_0)));
	gen double vio_1 = G*(1/(1+exp(-xbhat_1)));
	gen double lmv_0 = 1/(1+exp(-xbhat_lmv_0));
	gen double lmv_1 = 1/(1+exp(-xbhat_lmv_1));
	/* Generate predicted probabilities for democracy for each condition in both
	the split-sample and replication equations */	
	gen double vio_d1 = G*(1/(1+exp(-xbhat_d1)));
	gen double vio_d0 = G*(1/(1+exp(-xbhat_d0)));
	gen double lmv_d0 = 1/(1+exp(-xbhat_lmv_d0));
	gen double lmv_d1 = 1/(1+exp(-xbhat_lmv_d1));

	/* Plot the predicted probabilities of alliance violation across terroritorial 
	threat under varying conditions of change in leader's societal coalition and democracy */
	twoway  lowess vio_d1 terrthreat,  sort lpattern(solid) color(red) lwidth(thick)|| 
			lowess vio_d0 terrthreat, sort lpattern(solid) color(black)  ||
			lowess lmv_d1 terrthreat, sort lpattern(dash) color(red) lwidth(thick) ||
			lowess lmv_d0 terrthreat, sort lpattern(dash) color(black)		
					title("Democracy")
					ylabel(0(.01).04)
					xlabel(0(.05).42)
					scheme(s1mono) 
					ytitle("") xtitle("Territorial Threat") 
					legend(lab(1  "Democracy (w/ Relevance)") lab(2 "Non-Democracy (w/ Relevance)") 
							lab(3 "Democracy (Reduced)") lab(4 "Non-Democracy (Reduced)") rows(4) symx(*.5));
							
	graph save fig2a, replace;
	
	#delimit;
	twoway  lowess vio_1 terrthreat,  sort lpattern(solid) color(red) lwidth(thick)|| 
			lowess vio_0 terrthreat, sort lpattern(solid) color(black)  ||
			lowess lmv_1 terrthreat, sort lpattern(dash) color(red) lwidth(thick) ||
			lowess lmv_0 terrthreat, sort lpattern(dash) color(black) 
					title("{&Delta} in Leader's SC")
					ylabel(0(.01).04)
					xlabel(0(.05).42)
					scheme(s1mono) 
					ytitle("") xtitle("Territorial Threat") 
					legend(lab(1  "{&Delta} in Leader's SC (w/ Relevance)") lab(2 "No {&Delta} in Leader's SC (w/ Relevance)") 
							lab(3 "{&Delta} in Leader's SC (Reduced)") lab(4 "No {&Delta} in Leader's SC (Reduced)") rows(4) symx(*.5));
					
	graph save fig2b, replace;

	/* Graph combine predicted probability for Dem, Change in LSC and Territorial Threat*/
	#delimit;
	graph combine fig2a.gph fig2b.gph, scheme(s1mono) row(1) ycommon
					l1title("Pr(Violation)")
					note("Note: Estimates from Table 1, Models 2 and 3. All variables held at mean or median. Figure reports the" 
						 "predicted probabilities for the middle 95% of values of Territorial Threat from the samples.");
	graph save fig2, replace asis;


/* Beyond alliances: Trade application */	
/*Preamble*/
clear;
program drop _all;
clear matrix;
set seed 17;

	/* Calculate Rubin formula for uncertainty adding R-squared */
		capture program drop rubin;
		program define rubin;
			local names: colnames(Q1);
			/* MI param estimates */
			matrix Q = (Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q9 + Q10)/10;
			/* MI covariances     */
			matrix W = (W1 + W2 + W3 + W4 + W5 + W6 + W7 + W8 + W9 + W10)/10;
			/* MI R-squared */
			matrix r2 = (r21 + r22 + r23 + r24 + r25 + r26 + r27 + r28 + r29 + r210)/10;
			local  k = colsof(Q);
			forvalues i=1/10 {;
			  matrix QQ=Q`i'-Q;
			  if `i'==1 {;
				matrix B=(QQ)'*QQ;
				};
			  else matrix B = B + (QQ)'*QQ;
			};
			/* Covariance adjustment */
			matrix B=B/(10-1);		
			/* Final covariance   */
			matrix V=W+(1+1/10)*B;				
			ereturn post Q V ;
		end;
			
	/* Load trade data */		
	use NiemanGibler_demdiff_trade, clear;

	/* OLS for trade with democracy and territorial threat */
	forvalues i=1/10 {;
	areg ltrade BothDem_Threatmin`i' OneDem_Threatmin`i' BothDem OneDem threat_min`i' bothin onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island lareap comcol colony curcol comctry , a(year) robust cluster(dyad);
	matrix Q`i' = e(b);
		matrix W`i' = e(V);
		matrix  r2`i' = e(r2);
		nois _dots `i' 0;
	};

	rubin;
		estout , cells(b(star fmt(3)) se(par fmt(3)) ) starlevels(* .05) style(tex);
		matrix B = e(b);
		matrix V = e(V);
		matrix rename  r2 r2_threat;
		estimates store threat;
		
	/* Identify common observations using mean territorial threat score from the 10 draws */	
	areg ltrade BothDem_Threatmin OneDem_Threatmin BothDem OneDem threat_min bothin 
		onein gsp ldist lrgdp lrgdppc regional custrict comlang border landl island 
		lareap comcol colony curcol comctry , a(year) robust cluster(dyad);
	gen GN =1 if e(sample)==1;
				
	/* Figure 3 (and Figure A.4): Predicted Trade*/
	keep if GN==1;
	/* Collapse data to mean/median and expand for simulations */
	collapse (mean) threat=threat_min  (median) ldist lrgdp lrgdppc lareap bothin onein gsp  regional custrict comlang border landl island comcol colony curcol comctry (min) threat_min=threat_min (max) threat_max=threat_min;
	expand 100;	
		replace threat = threat_min + (_n/_N)*(threat_max - threat_min);
	expand 100;
	drawnorm b_BothDem_Threatmin b_OneDem_Threatmin b_BothDem b_OneDem b_threat_min b_bothin b_onein b_gsp b_ldist b_lrgdp b_lrgdppc b_regional b_custrict b_comlang b_border b_landl b_island b_lareap b_comcol b_colony b_curcol b_comctry b_cons, means(B) cov(V);
	/* Generate (unlogged) predicted values for two democracies, one democracy, and no democracy condititions */
	gen pr_trade_both = exp(b_BothDem_Threatmin*threat + b_BothDem + b_threat_min*threat + b_bothin*bothin + b_onein*onein + b_gsp*gsp + b_ldist*ldist + b_lrgdp*lrgdp +
		b_lrgdppc*lrgdppc + b_regional*regional + b_custrict*custrict + b_comlang*comlang +  b_border*border + b_landl*landl + b_island*island + b_lareap*lareap + 
		b_comcol*comcol + b_colony*colony + b_curcol*curcol + b_comctry*comctry + b_cons);
	gen pr_trade_one = exp(b_OneDem_Threatmin*threat + b_OneDem + b_threat_min*threat + b_bothin*bothin + b_onein*onein + b_gsp*gsp + b_ldist*ldist + b_lrgdp*lrgdp +
		b_lrgdppc*lrgdppc + b_regional*regional + b_custrict*custrict + b_comlang*comlang +  b_border*border + b_landl*landl + b_island*island + b_lareap*lareap + 
		b_comcol*comcol + b_colony*colony + b_curcol*curcol + b_comctry*comctry + b_cons);
	gen pr_trade_none = exp(b_threat_min*threat + b_bothin*bothin + b_onein*onein + b_gsp*gsp + b_ldist*ldist + b_lrgdp*lrgdp +
		b_lrgdppc*lrgdppc + b_regional*regional + b_custrict*custrict + b_comlang*comlang +  b_border*border + b_landl*landl + b_island*island + b_lareap*lareap + 
		b_comcol*comcol + b_colony*colony + b_curcol*curcol + b_comctry*comctry + b_cons);
	/* Calculate mean predicted values */	
	collapse (mean) pr_trade_both pr_trade_one pr_trade_none, by(threat);
	/*Plot mean predicted (unlogged) trade values for two, one, and no democracy conditions */	
		  twoway lowess pr_trade_both threat, scheme(s1color)
			lcolor(blue) 
			|| lowess pr_trade_one threat, 
			lcolor(red) lwidth(medthick) lpattern (longdash)
			|| 	lowess  pr_trade_none threat, 
			lcolor(black ) lwidth(medium ) lpattern(dash )
			ylabel(0(5000)15000)
			xtitle("Territorial Threat (weak link)")
			legend(lab(1  "Both Dem") lab(2 "One Dem") lab(3 "No Dem") symx(*.5) row(1));
			graph save fig3, replace asis;
	
log close;
